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Systematic reduction of sign errors in many-body problems: generalization of self-healing diffusion 

Monte Carlo to excited states 
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A recently developed self-healing diffusion Monte Carlo algorithm [PRB 79, 195 117] is extended to the cal- 
culation of excited states. The formalism is based on an excited-state fixed-node approximation and the mixed 
estimator of the excited-state probability density. The fixed-node ground state wave-functions of inequivalent 
nodal pockets are found simultaneously using a recursive approach. The decay of the wave-function into lower 
energy states is prevented using two methods: i) The projection of the improved trial-wave function into previ- 
ously calculated eigenstates is removed; and ii) the reference energy for each nodal pocket is adjusted in order to 
create a kink in the global fixed-node wave-function which, when locally smoothed, increases the volume of the 
higher energy pockets at the expense of the lower energy ones until the energies of every pocket become equal. 
This reference energy method is designed to find nodal structures that are local minima for arbitrary fluctuations 
of the nodes within a given nodal topology. It is demonstrated in a model system that the algorithm converges 
to many-body eigenstates in bosonic and fermionic cases. 

PACS numbers: 02.70.Ss,02.70.Tt 



I. INTRODUCTION 

Although several important chemical and physical proper- 
ties of matter are determined by the lowest energy electronic 
configuration (or ground state), a significant number of phys- 
ical properties are crucially dependent on the excitation spec- 
tra. These properties range from electronic optical excitations 
to transport and thermodynamic behavior 

While elegant theories that take advantage of the variational 
principle have been formulated for the ground state.fi3 the the- 
ories on the excitation spectra are far more complexi^ There- 
fore, although excited states are extremely important, our un- 
derstanding of them is limited as compared with the ground 
state. 

Diffusion quantum Monte Carlo (DMC) is the method of 
choice to obtain the ground state energy of systems with more 
than ^ 20 electrons. The DMC algorithm'* transforms the cal- 
culation of an excited state (e.g., the fermionic ground state) 
into a ground state calculation. The accuracy of the method 
depends, however, on a previous estimate of the zeros (nodes) 
of the wave-function. 

The ground state wave-function of most many-body Hamil- 
tonians 7i(R) is a bosonic (symmetric) wave-function with- 
out nodes. Any other eigenstate of a many-body Hamilto- 
nian 7Y(R) must have nodes in order to be orthogonal to the 
bosonic ground state. In the case of fermions (e.g., elec- 
trons), the ground state must be antisymmetric. Therefore, 
the electronic ground state is an excited state of the many- 
body Hamiltonian 7i(R) and must have nodes (hyper-surfaces 
in 3Ne space where the wave-function becomes zero and 
changes sign, being the number of particles). 

The standard diffusion Monte Carlo (DMC) approach^ 
finds the lowest energy E!^'^^'~^ of all the wave-functions 
that share the nodes 5t(R) of a trial wave-function 5'r(R), 
where R is a point in the 3A^e coordinate space. This low- 
est energy wave-function is denoted as the fixed-node ground 
state *^^(R). 



Since "no nodes" is a condition easy to satisfy, the ground 
state energy of a bosonic system can be found with a preci- 
sion limited only by statistical and time-step errors. For any 
other eigenstate 5'„(R), a good approximation of its nodal 
surface S'„(R) must be provided in order to avoid system- 
atic errors. Departures in S't(R) from the exact nodes 5„(R) 
cause, in general, errors of the energy as compared with the 
exact eigenstate energy.- For the fermionic ground state, the 
standard DMC algorithm provides only an upper bound of the 
ground state energy Moreover, if \E'„(R) is non degener- 
ate, any departure of St (R) from 5'„ (R) creates a kink in the 
fixed-node ground state.** Accordingly, accurate many-body 
calculations require methods to obtain and improve 5t(R). 
The problem of searching the exact nodes 5„(R), the sur- 
faces in 3A^e space where the wave-function of an arbitrary 
eigenstate n changes sign, is one of the outstanding problems 
in condensed matter theory. = 

This paper is the natural conclusion of earlier work. In 
Ref. MO we showed that even the exact Kohn-Sham- wave- 
functions cannot be expected to provide accurate nodal struc- 
tures for DMC calculations. However, we also showed that an 
optimal Kohn-Sham-like nodal potential exists. Subsequently 
in Ref. i8| we demonstrated that the nodes of the fermionic 
ground state wave-function can be found in an iterative pro- 
cess by locally smoothing the kinks of the fixed-node wave- 
function. We also showed that an effective nodal potential can 
be found to obtain a compact representation of an optimized 
trial wave-function with good nodes. While some details are 
rederived here, reading those papers before this one is highly^ 
recommended. 

In this paper the self-healing diffusion Monte Carlo method 
(SHDMC) is extended to find the nodes, wave-functions, and 
energies of low-energy eigen-states of bosonic and fermionic 
systems. 
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II. THE SIMPLE SHDMC ALGORITHM FOR THE 
GROUND STATE 

This paper describes how to extend the "simple SHDMC 
algorithm" (as described in Section III.C of Ref. [§1) to excited 
states. An extension to optimize the multi-determinant expan- 
sion, (see Section IV in Ref. ) is clearly possible and will be 
explained elsewhere. 

The ground state SHDMC algorithm builds upon the impor- 
tance sampling DMC method.^ The standard diffusion Monte 
Carlo approach is based on the Ceperley-Alder— equationt'^ 



df(R,T) 
dr 



V^/(R, r) - Vr (/(R, T)VR/n |vI.t(R) 

-[£;L(R)-£;T]/(R,r) , (1) 



where El{^) = [H*t(R)]/*t(R) is the "local energy," 
7i is the many-body Hamiltonian operator, R denotes a point 
in 3A^e space, and Et is a reference energy. Equation ^ is 
often solved numerically^ using a large number Nc of electron 
configurations (or walkers) which are points R; in the SiVg 
space. These walkers i) randomly diffuse according to the 
first term in Eq. ^ and ii) drift according to the second term 
a time 5t. In addition, iii) the walkers branch {or pass on} 
with probability p = 1 — exp[(£'L(R) — Et)5t] {or p = 
exp[(i?i(Ri) — i?T)(5r] — 1 }. To prevent large fluctuations in 
the population of walkers and excessive branching or killing, 
often a statistical weight is assigned to each walker A detailed 
review of the numerical methods used for minimizing errors 
and accelerating DMC calculations is given in Ref. [Tsl 

In the hmit of t ^ cx), the distribution function of the walk- 
ers in an importance sampling DMC algorithm is given byi 

/(R,r^cx)) = *5.(R)«';^Ar(R) e-^^?"""-^-)^ (2) 

N. 



lim lim ^^W/(j) j(R-R|) 



■I 

The R^ in Eq. (|2]i correspond to the positions of walker i at 
the step j for an equilibrated DMC run of Nc configurations. 
The original SHDMC method for the ground state was im- 
plemented in a mixed branching with weights scheme. For 
reasons that will be clear below, it is easier to formulate a 
method for excited states with a constant number of walkers 
with weights W.^ [k) which are given by 

W^'(fc) = e-NW-H'=*-, (3) 

with k being a number of steps, St the time step, and 



Ei{k)^lY.E,(R:>-'). 

e=o 



k-l 



(4) 



The energy reference Et in Eq. ^ is adjusted so that 
J2i (k) ~ assuming a constant Et for k steps. 

Note that setting all (fc) = 1 in Eq. ^ gives at equi- 
librium, by construction, a distribution /(R) = |^'t(R)P, 
because this is equivalent to setting El(R) = Et in Eq. ([T]i. 



If one sets the initial distribution of walkers as /(R, 0) = 
|5't(R)P, then the distribution of walkers at imaginary time 
T = kSr is given by 



/(R,r) = vl/y(R) 



-THf 



'«'t(R) 



(5) 



*t(R)*t(R, t) 



lim — y W/(fc)(5fR-RJ 



Therefore, at equilibrium and in a no branching approach, the 
weights VKf(fc) contain all the difference between /(R, r) 
and |*t(R)P . In Eq. ^ e^THpN jg tjjg fixed-node evo- 
lution operator, which is a function of the fixed-node Hamil- 
tonian operator 71 fn given by 

Ufn ='H-Et+oo lim 9{e- d„[S'T(R') - R]} . (6) 

The third term in the right-hand side of Eq. (|6|l adds an infinite 
potential at the points R with minimum distance to any point 
of the nodal surface (i,„[S'T(R') — R] smaller than eJ^ 
Using Eq. ^ one can formally obtain 



(R|*t(t)) =*T(R,r) = e 



-tHf 



*r(R) 



and using Eq. one obtains 

(RI^fat) = ^-i^wlR) = lim 1'T(R,T)e^ 



(7) 



(8) 



The trial wave-function ^[/^.(R) is commonly a product 
of an antisymmetric function $r(R) and a Jastrow'^ factor 
gj(R) Often $j'(R) is a truncated sum of Slater determi- 
nants orpfaffians $„(R): 



(9) 



In Ref. m we proved that we can evaluate e~'^^|5'T) for 
T ^ OO using a numerically stable algorithm. The analytical 
derivation of the algorithm*^ can be summarizedii here as 



= lim e 

T—>-QO 



lim nie-'^'^^e-^^^"")!*: 



f=0\ 
T / 



lim Y[{De-^'^'^ 



e=o\ 
T I 



(10) 
(11) 



T—^OO 



The operator D is defined in Eq. (fT6] l. Equation ( fTTT i means 
that the ground state |^'o)'^ can be obtained recursively by 
generating a new trial wave-function 1 4'^) from a fixed-node 
DMC calculation that uses the previous trial wave-function 
1 5'^"^), which is given by 



I* 



^ ^ ^ L> lim e-^^F« 



(12) 



FN I 
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Equation ( fT2] i means that new coefficients A„ of a truncated 
expansion of a trial wave-function of the form given in Eq. (|9]l 
are obtained numerically from the distribution of walkers of a 
DMC run as 



where 



(14) 



ani 



7(R) = 



-1 + ^1 + 2|v|2t 



with V 



(15) 



A complete explanation of our method is given in Ref. |8l 
Briefly here, our method systematically improves the nodes 
for three main reasons: 

1) The projectors in Eq. ( fT4l i include only functions $„(R) 
that retain all symmetries of the ground state. In more tech- 
nical terms, the ground state is expanded only with functions 
that belong to the same irreducible representation. This means 
that if the $„ (R) are determinants, for example, the bosonic 
ground state is excluded. Therefore, fluctuations that depart 
from the fermionic Hilbert space are filtered and do not prop- 
agate into the trial wave-function from one DMC run to the 
next SHDMC iteration. 

2) The projection of vI/j?^(R) into a finite set of $„(R) 
with low non-interacting energy can be shown^ to be equiv- 
alent to locally smoothing the kinks at the node of the fixed- 
node wave-function with a function of the form 



(R|^|R') = ~5 (R, R') = J2 '^nCRWni'R') 



(16) 



We proved that a large class of local smoothing functions have 
the same effect on the nodes as a Gaussian, under certain con- 
ditions, which includes the case of Eq. ( fT6] l. In turn, in Ref. 
m we proved that, to linear order in V&r', the convolution of 
a Gaussian with any continuous function has the same effect 
on the nodes as the imaginary time propagator e^^'^ [this 
allows replacing Eq. dTol i by Eq. (fTTt l. 

Thus our method can be viewed as the recursive application 
of two operators on the trial wave-function: i) g""^^^" that 
turns I^'t) into I^'fw) and ii) D that samples and truncates 
the expansion and changes the nodes as e^'^'^. Accordingly, 
our method is formally related to the shadow wave-function^i 
and the A-function approachJ^^i^ [see Eq. (fTotl. 

3) Finally, we argued that the method is robust against sta- 
tistical noise, because the kink should increase with the dis- 
tance between the exact node S'(R) and the node of the trial 
wave-function S't(R) [the kink must disappear for S't(R) = 
5(R)]. In addition, we took the relative error in A„ as trunca- 
tion criterion for D. 



III. EXTENSION OF THE SELF-HEALING DMC 
ALGORITHM TO EXCITED STATES 

A detailed explanation of the advantages and limitations 
of the standard fixed-node approximation for excited states is 
given in Ref.jsjThis paper explores the possibility of overcom- 
ing these limitations in calculating excited states by excluding 
the projection of lower energy states from the set of ^„(R). 
However, in to follow this path the problem of inequivalent 
nodal pockets has to be addressed. 



A. Inequivalent nodal pockets 

The expression "nodal pocket" denotes a volume in 3Ne 
space enclosed by the nodal surface S't(R). It has been 
shown^ that the ground state of any fermionic Hamiltonian 
with a local potential has nodal pockets that belong to the 
same class, meaning that the complete 3Ne space can be cov- 
ered by applying all symmetry operations (e.g., particle per- 
mutations) to just one nodal pocket. Therefore, if the trial 
wave-function is obtained from such a Hamiltonian, all nodal 
pockets are equivalent by symmetry. For the ground state, one 
can obtain the fixed-node wave-function in just one pocket and 
map it to the rest of the 3Ne space using permutations of the 
particles and other symmetries of H. 

In the case of arbitrary excited states, there are inequiva- 
lent nodal pockets that present a challenge to the fixed-node 
approach;=2 Due to this inequivalent pocket problem, alter- 
natives to the fixed-node method and variations have been 
^■ied,2i.22.23,24,25,26,27.28.29 Self-heaUng DMG** implicitly takes 

advantage of the equivalence of nodal pockets in the fermionic 
ground state and must be extended to the inequivalent pocket 
case. For this reason a nonbranching formulation is used in 
the excited state case. 



B. Equilibration of walkers in inequivalent nodal pockets 

A first complication, which has a simple solution, of the 
nonbranching fixed-node approximation is that the number of 
walkers in each nodal pocket is also fixed by the nodes. As 
a consequence of the drift or "quantum force" term [second 
term in Eq. ([T])], the walkers are repelled from the regions 
where the wave-function is zero and they cannot cross the 
node for 6t 0. The fact that the population in each nodal 
pocket is fixed has no consequence for the ground state be- 
cause all nodal pockets are equivalent. For the ground state it 
is not important in which nodal pocket the walker is trapped 
because particle permutations can move every walker into the 
same nodal pocket and the projectors ^„ (R) in Eq. ( [T4l i are 
invariant under such permutations. 

However, in the case of excited states, which have more 
nodes than those required by symmetry^^ there are inequiv- 
alent nodal pockets. In a nonbranching DMC scheme with 
weights, the population is locked from the start in a set of 
pockets. If the initial distribution of Nc walkers is chosen 
with a Metropolis algorithm to match |5't(R)P, there would 
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be random variations in the starting population of the order 
of \J Nc/Np, where Np is the number of inequivalent nodal 
pockets. This would cause systematic errors if the wave- 
function coefficients A„ were sampled without taking preven- 
tive measures. Moreover, even if the initial numbers of walk- 
ers in each pocket were set "by hand" (to be proportional to 
the integral |'I't(R)P in each pocket), the resolution of the 
sampling cannot be better than \/Nc- The importance of this 
error grows if Nc is small or if the number of inequivalent 
nodal pockets is large. 

To prevent this error from occurring, some walkers are sim- 
ply allowed to cross the node after the wave-function coef- 
ficients are sampled. At the end of a sub-block of k steps, 
for every walker i at R^, a random move AR^ is gener- 
ated with a Gaussian distribution using = St', without 
the drift velocity contribution. This move is accepted only 
if the wave function changes sign with a Metropolis prob- 
ability p = max{l, [*T(Ri+AR)/*T(Rj)]^}- This en- 
sures that i) the distribution of walkers remains proportional 
to |\1/t(R)P and ii) the average number of walkers in each 
pocket is proportional to the integral of |'I'7'(R)p as the num- 
ber of sub-blocks M tends to oo. 



D. Direct projection 

While the trial wave-function can be forced to be orthogo- 
nal to the ground state, or any other excited state calculated 
before, the fixed-node wave-function can develop a projection 
into lower energy states, because the DMC algorithm only re- 
quires \I'i?Ar(R) to be zero at the nodes 5t(R). To prevent 
excited states from drifting into lower energy states, let me as- 
sume, for a moment, that approximated expressions of the ex- 
cited states (R|e-^|$„) = *„(R) = e"'(^)$„(R) withn < ly 
can be obtained and used to build the projector 



P 



1-^|$„)(C| 



(17) 



where the operator e"^ is the multiplication by a Jastrow. Since 
the |$„) shall be obtained statistically, they will have errors 
and will not form an orthogonal basis in general. There- 
fore, ($*| are elements of the conjugated basis that satisfy 
(•S* I'I'm) = Sn.m- They Can be constructed inverting the Over- 
lap matrix Sn,,n = (^'ril^'m) aS 



C. Unequal fixed-node energies in inequivalent nodal pockets 



(18) 



A second complication of the fixed-node approach for the 
general case of excited states appears because small depar- 
tures of 5t(R) from the exact nodes S'„(R) often will re- 
sult in inequivalent nodal pockets having fixed-node solutions 
with different fixed-node energies. When nodal pockets are 
not equivalent, a standard DMC algorithm will converge to a 
"single nodal pocket" population. In this case, the lowest en- 
ergy pocket will contain all the walkers in a branching algo- 
rithm [or all significant weights (VF/(fc) ^ )]. Accordingly, 
the average energy sampled will correspond to the lowest en- 
ergy nodal pocket, which will be different from that of the true 
excited-state energy (see Chapter 6 in Ref. and references 
therein). 

If the coefficients of an excited-state fixed-node wave- 
function are sampled with the same procedure used for the 
ground stated [see Eq. (fT3]l1. they would correspond to a func- 
tion that is different from zero just at the class of nodal pock- 
ets with lowest DMC energy and zero everywhere else. This 
function will not be, in general, orthogonal to the lower en- 
ergy states. Moreover, this will result in kinks at the nodes 
in the wave-function sampled with Eq. ( fT3] l between lowest 
energy nodal pockets and inequivalent ones. 

A first preventive measure to avoid a single pocket popula- 
tion is to avoid the limit r ^ oo in Eqs. ( fTTT i and ( fT2] i which 
replaces l^p^) by c^'^^^'^fn in Eq. ([T2]i. As a re- 

sult k in Eq. ( fT3T l is limited to small values, which brings all 
values of closer to 1. Since the approach is recursive, 

the limit of r ^ oo is reached as £ ^ oo (since successive 
applications of the algorithm are accumulated in |5'f-)). In 
addition, to prevent the wave-function from falling into lower 
energy states, two techniques are used: i) direct projection and 
ii) unequal reference energies. 



Then, the extension of the self-healing algorithm to the next 
excited can be rederived analytically as follows: 



I* 



lim P e~^^P|«'^=° i) 

r — >OQ 



lim P J] L-^^r'+k5r)np\ |^^=0^^ 



lim P 



lim P 



n(. 

t 

no 



-kSrH 



0(1-1) 

FN P 



P 



T,i/+1/ 



(19) 



Equation ( fT9] l means that for any initial trial wave-function 
I*t"°+i) with P\^^T^^^) ^ 0, one can obtain the next ex- 
cited state l^'^+i) recursively. The numerical implementation 
of the algorithm for excited states (see Section|lV]for details) 
is almost identical to the ground state version** with three dif- 
ferences: i) there is no branching and the product kSr is cho- 
sen so as Wl{k) ~ 1 [see Eq. (fT3])1, ii) the projection of 
the vector of coefficients A„ into the ones corresponding to 
eigenstates calculated earlier is removed with P, and iii) some 
walkers cross the node after k time steps (see above). 

Eq. (Hg holds in the limit of oo, St ^ 0, St' 0, 

£kST oo, and £St' ^ oo. In the derivation of Eq. iT% . the 



following properties were used: P^ = P, and [H, P] 
Ref. d it was shown that, under certain conditions. 



0. In 



S 



S 



pm 



(20) 
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that is, the nodes of the two functions in the brackets are ap- 
proximately the same. 

Note that the second term in brackets of Eq. (fTTl i has pre- 
cisely the form given in Eq. ( fT6b . By construction, this term 
would generate a function with nodes corresponding to a lin- 
ear combination of lower energy eigenstates. The projector P, 
instead, excludes any change in the wave-functions introduced 

by the projection and sampling operator D or by e '^"■pn in 
the direction of lower energy wave-functions (which includes 
their nodes). 



E. Adjusting the reference energy in each nodal pocltet 

If walkers at one side of the node have more weight than at 
the other (because of inequivalent pockets with different fixed- 
node energies), the propagated wave-function obtained by 
sampling the walkers will be multiplied by a larger (smaller) 
factor for the low (high) energy side of the nodal surface. This 
generates an additional contribution to the kink at the node 
that, when locally smoothed, increases the volume of lower 
energy pockets at the expense of the higher energy ones, caus- 
ing the volume of the lower (higher) energy pockets to grow 
(diminish). This, in turn, will have an impact on the kinetic 
energy: due to quantum confinement effects, the difference in 
fixed-node energies will increase in the next iteration. This 
very interesting effect in fact acts to our advantage by helping 
us to find the ground state even when starting from a very poor 
wave-function.— For excited states, this effect is prevented by 
i) limiting the maximum value of k and ii) the projector P in 
Eq. dTgl l. However, the eigenstates will have statistical 
errors that can create systematic errors in the higher states. To 
partially prevent these errors, and to limit the number of or- 
thogonality constraints, the energy reference can be changed 
in order to invert this contribution to the kink to our advantage. 

While a single reference energy Et can still be used for 
the DMC run in each block, the projectors of Eq. ( fT3] l are 
redefined using a reference energy dependent on the nodal 
pocket. In addition, following a suggestion of C. Umrigar,^' 
the change in the coefficients 5A„ is sampled instead of the 
total value A„ . 



the nodal pocket where the walker i is trapped; however, only 
the last two-thirds of the block are used to accumulate values 
to allow (jo) to equilibrate. 

It was argued before that, for (3 = 0, the differences in 
the fixed-node energies of neighboring nodal pockets create 
a contribution to the kink that, when locally smoothed, in- 
creases the volume of nodal pockets with low fixed-node en- 
ergy. For (3 > 1, it is likely that this contribution to the 
kink is inverted so that the volume of the lower (higher) en- 
ergy pockets is reduced (increased) by the smoothing func- 
tion ( fT6l ). Therefore, it can be assumed that a value of /3 > 1 
should stabilize the higher energy nodal pockets, increasing 
their volume and, thus, reducing their energy. This process 
will stop when the fixed-node energy of all nodal pockets be- 
comes equal. 

Note that by introducing this artificial contribution to the 
kink, one may stabilize some nodal structures, preventing 
nodal fluctuations that reduce the energy of one nodal pocket 
at the expense of the others. However, fluctuations that lower 
the energy of every nodal pocket are not prevented. Therefore, 
if several eigenstates have the same nodal topology, higher en- 
ergy states could drift into lower energy ones if orthogonality 
constraints [see Eq. ( fTTt l are not imposed. 

Finally, note that choosing /? > 1 can also cause problems 
if the quality of the wave-function is not good or if the statis- 
tics is poor For example, a small statistical fluctuation in the 
values of A„ could create a new nodal pocket with high en- 
ergy. In successive blocks (as £ increases), this pocket will 
grow at the expense of the others, causing the total energy to 
rise. 



IV. SHDMC ALGORITHM FOR EXCITED STATES 



A„ — A 



— T 

2—1 



{SXn) 



where (3 is an adjustable parameter and 



A basis of <&„ (R) must be constructed, taking advantage of 
all the symmetries of 7i 3 The $„ (R) should be selected to 
be eigen-functions of a noninteracting many -body system^, be- 
longing to the same irreducible representation for every sym- 
metry group of H. The calculation must be repeated for each 
irreducible representation. Note that the same algorithm is 
used for bosons or fermions: the only difference is the basis 
(21^d to expand the wave-functions. 

The calculation of excited states with SHDMC is composed 
1) £,* (R"^ ) 7(R''?^ ^ sequence of blocks. Each block i has M sub-blocks with 
" * K standard DMC steps. 

The basic algorithm is the following: 



(22) 



is the weighted average of the local energy during the lifetime 
of the walker i since the start of the block or the last time it 
crossed the node at step Jq. If (3 = 1 is selected in Eq. ( |2TI ). 
the factor e~'^^^'^~^i ^^"^1 just replaces in the definition of the 
weights [see Eq. ©J Et by Ej {Jq). The energy Ej (jo) for 
j^jo ^ k is expected to converge to the fixed-node energy of 



1. An initial set of coefficients for the expansion of the 
trial wave-function is selected. 

2. The changes (5A„ are accumulated [see Eqs. ( fT4l ) and 
( |2TI )1 at the end of each sub-block. Some walkers near 
the node can cross it at the end of each sub-block. 

3. At the end of each block £, the error in dXn is evaluated. 
If this error is larger than 25% of A„ + ^A„, then A„ is 
set to zerOfS otherwise, A„ is set to A„ + SXn- 
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4. A new trial wave-function is constructed at the end 
of each block I using the new values of the coeffi- 
cients sampled after removing with P the projection 
into eigenstates calculated earlier. 

5 . If the scalar product between the vector of new 5 A„ with 
the one obtained in the previous block — 1) is posi- 
tive, the number of sub-blocks M is increased by one. 
Otherwise, M is multiplied by a factor larger than one 
(e.g., 1.25). This factor increases the statistics reducing 
the impact of noise 

6. Steps 2-6 are repeated until the variance of the weights 
Wlik) is smaller than a prescribed tolerance (see Fig. 
|6]in SectionlV). 

7. The projector P is updated to include the new excited 
state. 

8. Steps 1-7 are repeated until a desired number of excited 
states is obtained. 



A. Remarks 

Some points about the application of the algorithm should 
be addressed before discussing the results. 

• In this paper, to test the method, intentionally poor trial 
wave-functions have been selected as a starting point. 
Good initial wave-functions and a good Jastrow are ad- 
vised in real production runs in large systems. Methods 
to select good initial trial wave-functions will be dis- 
cussed elsewhere. 

• Time-step errors and, in particular, persistent walker 
configurations^^ can cause significant problems. When 
this happens it often results in an increase in the er- 
ror bar of every A„ which causes a large reduction in 
the number of coefficients retained in the trial wave- 
function. This problem is avoided in the algorithm by 
discarding the entire block if a 50% reduction in the 
number of basis functions retained is detected. Never- 
theless, if the quality of the initial (R) is bad, it is 
strongly recommended to reduce the time step Jr. As 
the quality of the wave-function improves with succes- 
sive iterations, one can increase 8t. For fast conver- 
gence Vk&T should be of the order of the interparticle 
distance. 

• As a strategy, it is better to run at first using /? = in 
Eq. ( I2TI 1 including every state calculated before in P 
[see Eq. (fTTIil. Once the wave-function \I't(R) is con- 
verged, one can set P = 1 and [3=1 and monitor if 
^frCR) evolves into a subset of lower energy states. To 
prevent the propagation of errors of every lower energy 
state included in P into the next excited state, a run in- 
cluding only this subset in P can be performed. 

• To obtain accurate total energies, a long run with large 
k is required (this is almost a standard DMC run). 



• SHDMC should not be used blindly as a library routine. 
The calculation of excited states with SHDMC is a task 
that will probably remain limited to quantum Monte 
Carlo experts. While, in contrast, density functional ap- 
proximated methods have suddenly become very easy 
to use, it is not quite clear to the author that requir- 
ing expertise and a deep understanding is a disadvan- 
tage. Any new code using SHDMC should be tested 
in a small system where analytical solutions or results 
with an alternative approaches, are available. The com- 
parison with a soluble model is presented in the next 
section. 



V. APPLICATIONS TO MODEL SYSTEMS 

This section compares the methods described above for the 
calculation of excited states with SHDMC, with full configu- 
ration interaction (CI) calculations in the model system used 
in Refs. andi. 

Briefly, the lower energy eigenstates are found for two elec- 
trons moving in a two dimensional square with a side length 1 
with a repulsive interaction potential of the forroi^ V{r, r') = 
87r^7Cos [Q!7r(a:; — x')] cos [a7r(i/ — ?/')] with a = I/tt and 
7 = 4. The many-body wave-function is expanded in func- 
tions $„ (R) that are eigenstates of the noninteracting system. 
The $„(R) are linear combinations of functions of the form 
sin{in^TTx,^) with < 7. Full CI calculations are per- 
formed to obtain a nearly exact expression of the lower energy 
states of the system *„(R) = Y.m am*™(R)- 

We solve the problem both for the singlet and the triplet 
case. The singlet state of this system is bosonic-like, since the 
ground state wave-function has no nodes. The lowest energy 
excitations of the noninteracting problem $„(R) that have 
the same symmetry (that is, that are invariant under exchange 
of particles, and under all symmetry operations of the group 
D4) are selected to expand H. For the case of the triplet, the 
wave-function must change sign for permutations of the par- 
ticles. The ground state is, however, degenerate (belongs to 
the E representation of D4). The E representation can be de- 
scribed by wave-function even (odd) for reflections in x and 
odd (even) for reflections in y. We choose the wave-functions 
that are odd in the x direction: belonging to a D2 subgroup of 
the D4 symmetry. For more details on the triplet ground state 
calculations, see Refs. [lol and [si 

To facilitate the comparison with the full CI results, projec- 
tors (R) are constructed with the same basis functions used 
in the CI expansion. For the same reason, no Jastrow function 
is used [ J = in Eq. (fT4]i1. 

To test the method, poor initial trial wave-functions are in- 
tentionally chosen as follows: For the ground state the lowest 
energy function of the noninteracting system is selected. For 
the 71*'* (n = iy+1) excited state, the initial trial wave-function 
I'^Tn') constructed by completing the first v columns of 
a determinant with the first ly + 1 coefficients of the 1/ eigen- 
states calculated before. Subsequently, the vector of cofactors 
of the last column was calculated. The coefficients of this vec- 
tor are used to construct a trial wave-function orthogonal to all 
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FIG. 1 ; (Color online) Self-healed DMC run obtained for successive 
eigenstates belonging to the Ai (trivial) irreducible representation of 
the group D4, in the singlet state. Black lines denote the average 
value of the local energy. The horizontal blue dashed lines mark the 
energy of the corresponding excitation in the full CI calculation. 



the eigenstates calculated earlier. 

Figure [T] shows the results of successive SHDMC runs for 
the singlet ground state and the next 8 excitations that belong 
to the same symmetry (total spin S ~ 0, and irreducible rep- 
resentation Ai of the group D4). The SHDMC calculations 
were done using Nc = 200 walkers with a sub-block length 
fc = 50, a time step St = 0.0002^ 6t' = 0.002 (for the 
ground state St' = ) and, /3 = 1 in Eq. i2T[ . 

The lines in Fig. [T]join the values obtained for the weighted 
average of the local energy El{II) for each time step. The 
horizontal dashed lines mark the energy of the nearly ana- 
lytical result obtained with full CI. The agreement between 
SHDMC and full CI is extremely good. As higher energy 
eigenstates are calculated however, and the number of nodal 
pockets and nodal surfaces increases, time step errors start to 
play a dominant role. In particular, for the 9*^ excitation (not 
shown) St must be reduced. 

The occasional peaks (or drops) observable in the data are 
correlated with the update of \I't(R), and their reduction also 
reflects a systematic improvement in the trial wave-function. 
At the end of each block, the trial wave-function coefficients 
Xn are updated and all weights are reset to 1. They gradu- 
ally reach equilibrium values when new energies are sampled, 
completing a sub-block of length k. As a result, at the be- 
ginning of each block, the energy sampled is the average of 
the trial wave-function energy, which is often different than 
the DMC energy sampled thereafter (but it can be smaller or 
higher for a bad trial wave-function with small Nc). 

One interesting result is that some orthogonality constraints 
are not required to obtain some excited states. This is the case, 
for example, of the first excited state calculated with /3 = 1. 
This is presumably due to the fact that the number of nodal 
pockets is different for the excited state and the ground state 
and the decay path from the first excited state to the ground 
state is obstructed by the formation of a kink between inequiv- 



TABLE I: Values obtained for Lrp [see Eq. l |23t ] for a total of (a) 4 x 
lO"* (b) 8x10'' and (c) 12 x lO"* DMC steps and corresponding eigen- 
energies for two electrons in a square box with a model interaction. 
The logarithm of the residual projection Lrp of the SHDMC wave- 
function with the corresponding full result CI is given for different 
eigenstates belonging to the same symmetry of the ground state as 
a function of the number of steps used to sample the wave-function. 
The states are included in the order they were obtained. 



State 


Spin 


Rep. 


Lrp 


Lr-p 


Lrp 


CI 


SHDMC 










a 


b 


C 


Energy 


Energy 







S 


Ai 


-14.84 


-15.05 




328.088 


328.089 


(2) 


1 


S 


Ai 


-6.80 


-8.85 




374.106 


374.103 


(6) 


2 


S 


Ai 


-7.23 


-8.69 




409.960 


409.954 


(3) 


3 


S 


Ai 


-4.42 


-6.07 




418.508 


418.66 


(2) 


4 


S 


Ai 


-3.65 


-5.01 




454.630 


454.84 


(2) 


6 


S 


Ai 




-4.85 


-6.22 


477.019 


477.100 


(5) 


7 


S 


Ai 


-3.90 


-5.26 




492.216 


491.98 


(1) 


5 


S 


Ai 


-5.60 


-6.17 




468.854 


468.845 


(13) 


8 


S 


Ai 


-5.09 


-6.49 




503.805 


503.92 


(1) 





T 


E 


-8.49 


-8.71 




342.137 


342.191 


(5) 


1 


T 


E 


-4.37 


-4.35 




385.908 


387.80 


(1) 


3 


T 


E 


-3.06 


-3.35 




422.670 


423.60 


(2) 


5 


T 


E 


-4.04 


-5.48 




438.791 


438.70 


(1) 


2 


T 


E 


-2.31 


-2.31 




411.887 


416.07 


(1) 



alent nodal pockets if a value of /? « 1 is used. This is also 
the case for states 6 and 7, which were obtained before state 5 
despite the fact that they have higher energy. 

A similar effect is observed in some triplet excitations. Due 
to the choice of initial trial wave-function and the kink in- 
duced by /3 = 1, the 3'"'* excitation is found before the 2"'*, 
and the 5*'* is obtained before the 2"'' and the 4*'' . This inter- 
esting effect disappears if /3 = is chosen. 

TableUshows the logarithm of the residual projection 

L,p = log(l-|(*,'^^|*„)|) (23) 

of the excited state wave-function |4'„) sampled with 
SHDMC onto the corresponding full CI result as a 

function of the number of iterations for different eigenstates. 
The states are ordered as they first appear in the calculation. 

In addition. Table U compares the values of the eigen- 
energies obtained with CI and SHDMC. The agreement is 
very good. In some cases the difference is larger than the er- 
ror bar. This might signal that small nodal errors remain. Note 
that there is no upper bound theorem for excited states but for 
the ground state within an abelian irreducible representation^ 

Figure|2]shows L,.p at the end of each block for the ground 
state and low-lying excitations of the system as a function of 
the total number of SHDMC steps. The calculations were 
done by first running ^ 40 000 SHDMC steps for each eigen- 
state before starting the calculation of the next. Subsequently, 
an additional set of ~ 40 000 SHDMC steps was run, improv- 
ing the projector P. The kinks in the data around ~ 40 000 
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FIG. 2: (Color Online) Logarithm of the residual projection [see Eq. 
l |23H for the ground (square), first (diamond), second (up triangle) 
and third (down triangle) eigenstates with Ai symmetry and S=0. 



are due to the changes in the coefficients of the lower energy 
states involved in P [see Eq. (fTTIil. 

One important conclusion of Table I] and Figure |2] is that 
errors in the determination of lower energy states calculated 
earlier only propagate "locally" because of the orthogonality 
constraints in Eq. (fTTj i. This error does not have a strong im- 
pact on much higher energy excitations. This is apparently 
due to the fact that each newly calculated excitation tends to 
occupy the Hilbert space left by lower excitations due to sta- 
tistical error. This is clear, for example, for the 5*'* and 8*'' 
excitations, which have an error much smaller than several 
excitations calculated earlier (e.g., 3'"'* and 4*''). The error in 
the 3'"'' and 4*'' excitations is mainly due to mixing among 
themselves. This result is important because it means that the 
present method can be used to calculate several higher excita- 
tions in spite of the errors in lower energy ones. 

Figure [3] shows the evolution of the values of the coeffi- 
cients Afj of l^*^) as a function of the coefficient index n for 
the 5*'* excited state corresponding to the singlet configuration 
of the Ai representation of the group D4. The shade of gray is 
light for the older (small £) coefficients and deepens to black 
for the final results (large £). The calculation started from a 
trial wave-function orthogonal to the states calculated before 
as described above. 

The coefficients of the wave-function sampled with 
SHDMC overlap with the ones obtained with full CI (see Ta- 
ble nil. Similar results are obtained for all the other excited 
states calculated. An important observation is that the co- 
efficients A„ evolve continuously towards the exact solution, 
which suggests the possibility of accelerated algorithms that 
extrapolate the values of (5A„. 

Some eigenstates are significantly more difficult to calcu- 
late than others. This is typically the case for eigenstates 
with similar eigenvalues (e.g., the 6*'' excitation in the sin- 
glet case). A bigger challenge, however, is when £'l(R) is ill 
behaved, for example, the case of the 2"'', 4*'', and 6*'' excita- 
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5'^ excitation (singlet A,) 

SHDIVIC 
♦ Full CI 
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Wave-Function coefficient index n 



FIG. 3: (Color online) Change in the values of the multi-determinant 
expansion as the DMC self-healing algorithm progresses for the 5"* 
excited state of the singlet state of j4i symmetry. Light gray col- 
ors denote older coefficients, whereas darker ones denote more con- 
verged results. The full CI results are highlighted in small red dia- 
monds. 



tions of the triplet state. Even the full CI wave-function with 
300 basis functions has a large variance for ElCR). In that 
case the coefficients obtained with SHDMC and CI are differ- 
ent. This is due to the fact that the two methods minimize dif- 
ferent things; CI minimizes ( — £'„)^|\I'„) on a trun- 
cated basis, and SHDMC minimizes J EL(R) f{R,T)dR 

with {^t\P\^t) = (*r|*r)- Accordingly, the fact that the 
results are different indicates that neither calculation, CI or 
SHDMC, is converged with the basis chosen. The 4*'' and e*'' 
excitations with E symmetry in the triplet case obtained with 
SHDMC are a linear combination of the corresponding ones 
in full CI. 

Figure |4] shows the effect of P and (3 [see Eq. d^TT il on a 
SHDMC run. The figure shows the average of the local en- 
ergy El (R) for two calculations that start from the final trial 
wave-function obtained for the 8*'' singlet excitation with Ai 
symmetry (please compare it with Fig. [T]!. Both calculations 
were run with the same parameters as in Fig. [T]with two ex- 
ceptions: i) P = 1 was used, which removes the orthogonality 
constraints, and ii) one calculation was run with (3 = 1.05 and 
the other with /3 = in Eq. (I2TI 1. An initial number of blocks 
M = 20 was used. 

Both calculations depart from the initial configuration. 
However, the run with /3 = falls very quickly to the sin- 
glet ground state. The calculation with (3 = 1.05 remains 
much longer in the vicinity of the 8*'* excitation. This clearly 
shows the stabilizing effect unequal energy references on ex- 
cited states. Since presumably the 8*'* excitation is not the 
minimum of its nodal topology, it finally drifts away. For the 
(3 = 1.05 case with 5t = 0.0002, the algorithm becomes nu- 
merically unstable to noise after the ^ 50] , 000 time step be- 
cause the variance in the distribution of weights of the walkers 
increases and the statistics is dominated by a reduced number 
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FIG. 4: (Color online) Average of the local energy El (R) as a func- 
tion to the number of DMC time steps for two SHDMC runs with 
P = 1 starting from a converged trial wave-function corresponding 
to the 8"' singlet excitation of Ai symmetry with a) 13 — 1.05 and 
b) /? = in Eq. \2U . The dotted lines mark the beginning of some 
of the fixed-node DMC blocks of a SHDMC run for the /3 = case. 
Same conventions as in Fig. [T] 
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FIG. 5: (Color online) Average of the local energy El(R.) of 200 
walkers as the SHDMC algorithm converges to the ground, first and 
second eigenstates with Ai symmetry and S=0 of two electrons with 
Coulomb interactions in a square box. 



of walkers. 

In contrast, the first excitation does not drift with /3 ~ 1 and 
P = 1 (not shown). 



A. Coulomb interaction results and discussion 

The use of a simplified electron-electron interaction facili- 
tates the CI calculations and the validation of the optimization 
method. However, it is also important to test the convergence 
and stability of the method with a realistic Coulomb interac- 
tion as in the case of the ground statCf^ 



10 15 20 

SHDMC block index (^) 



FIG. 6: (Color online) Logarithm of the variance of the weights of 
the walkers distribution as a function of the SHDMC block index £ 
for the 2"'' excitation with A^ symmetry with Coulomb interaction 
(see Fig. Is}. The lines are visual guides. 



The results shown in this section have an interaction poten- 
tial of the formic V{r, r') = 207rV|r - r/| as in Ref. .8. To 
mimic the difficulties that the algorithm would have to over- 
come in larger or more realistic systems, the Jastrow term 
is not included, i.e. J = 0. Most SHDMC parameters are 
the same as in the model interaction case. All calculations 
with Coulomb interactions were run with [3^0, the initial 
number of sub-blocks AI = 6, and the time step reduced to 
6t = 0.0001. The initial trial wave-functions were selected 
with the criteria used for the model case. 

Figure |5] shows the average of the local energy El(R) ob- 
tained for the ground state and the first two excitations with 
the same symmetry (singlet Ai). The results are qualitatively 
similar to those obtained with the model potential. It is evi- 
dent from the data that the variance of El (R) and its average 
are reduced as the wave-function is optimized. Occasionally, 
El(R.) might rise when P is updated (improving the descrip- 
tion of lower energy states). 

The energy of the singlet ground state is 400.749 ± 0.013, 
which is only slightly smaller than the lowest triplet energy^ 
402.718 ± 0.008 with symmetry E. These energies are very 
close because of the dominance of the Coulomb repulsion as 
compared to the kinetic energy, which forces the particles to 
be well separated and therefore the cost of a node in the triplet 
state is small. This result is consistent with the choice of pa- 
rameters that sets the system in the highly correlated regime. 
The energies obtained for the first and second excitations ar»i^ 
468.56 ± 0.09 and 515.50 ± 0.08 respectively. 

While Figs. [T] and |5] are qualitatively similar, the results 
shown in Fig. [T]are more convincing since they are directly 
compared with full CI calculations and they are less noisy, as 
noted by one referee. When the model interaction potential is 
replaced by a Coulomb interaction, full CI calculations are 
still possible, but they involve the numerical calculation of 
16471 integrals with Coulomb singularities. CI calculations 
are typically done using a Gaussian basiSj-'-^ which limits the 
impact of the matrix element integrals of these singularities. 
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However, as the size of the system increases, CI calculations 
become too expansible numerically. Accordingly, self-reliant 
methods to validate the quality of the SHDMC wave-functions 
must be developed. 

As noted earlier, in a fixed population scheme, the weights 
contain all the difference between /(R, t) and |'J't(R)P ■ 
Since /(R, r) and |vI/j.(R)|2 should be equal if 'I't(R) is an 
eigenstate, the variance of the weights can be used to measure 
the quality of the wave-function. Figure|6]shows the evolution 
of the logarithm of the variance L^ar of the weights of the 
walkers (k) [see Eq. ©] as a function of the SHDMC 
block index £. L^ar is evaluated as 



i.a. = log^— (24) 

By using a linear order expansion in St in Eq. (O and using 
Eq. (Ell, it is straightforward to relate Eq. (l24b to the variance 
of Ei'- (fc). The latter is an average of ^^^(R). A common 
measure of the quality of the ground state wave-function is 
the variance of (R). 

The results shown in Fig. |6] correspond to the 2'"^ singlet 
excitation with Ai symmetry (see Fig. |5]). Similar results 
are obtained for the ground state and the first excitation (not 
shown). The error bar in Lyar is smaller than the size of the 
symbols. The fluctuations in Lyar result from the random 
fluctuations of the coefficients A„ that are obtained statisti- 
cally. Note that in spite of the noise, a clear trend shows the 
improvement of the quality of the wave-function and Et as 
the SHDMC algorithm progresses. However, these improve- 
ments are not uniform, which is reflected by the oscillations 
in Lyar in Fig. |6]and in the amplitude of ^^^(R) in Fig. |5] 
A careful user of SHDMC should track L^ar and use the best 
quality wave-function to calculate energies and P. 



VI. SUMMARY 

An algorithm to obtain the approximate nodes, wave- 
functions, and energies of arbitrary low-energy eigenstates 
of many-body Hamiltonians has been presented. This algo- 
rithm is a generalization of the "simple" self-healing diffu- 
sion Monte Carlo method developed for the calculation of the 
ground state of fermionic systems^^ which in turn is built upon 
the standard DMC method^ 

At least in the case of the tested system, wave-functions and 
energies that continuously approach fully converged configu- 
ration interaction calculations can be obtained depending only 
on the computational time. The wave-function, in turn, allows 
the calculation of any observable. 

It is found that some special eigenstates, presumably the 
minimum energy eigenstate for a given nodal topology, can 
be obtained without calculating the lower excitations by ar- 
tificially generating a kink in the propagated function using 
unequal energy references in different nodal pockets. 

The present method can be implemented easily in existing 
codes. Ongoing tests on the ground state method^ in larger 



systems give serious hopa** that the current generalization 
will also be useful. 

While there are methods to obtain the excitation spectra of 
a many-body Hamiltonian in a variational Monte Carlo con- 
text^**'^^^ they require obtaining the Hamiltonian and the over- 
lap matrix elements. This requirement would present a chal- 
lenge for very large systems. SHDMC is a complementary 
technique that could potentially scale better for larger sizes. 
The evaluation and storage of the matrix elements of H is not 
required. The number of quantities sampled [the projectors 
^„(R), Eq. ( fT4l i1 is equal to the number of basis functions ??{,. 
In contrast, energy minimization methods or configuration in- 
teraction (CI) require the evaluation of matrix elements. 
In addition, the solution of a generalized eigenvalue problem 
with statistical noise is avoided. This can be an advantage in 
very large systems since algorithms for eigenvalue problems 
are difficult to scale to take maximum advantage of large su- 
percomputers. In contrast, the sampling of a large number 
of determinants can be trivially distributed on different pro- 
cessors. Moreover, recent advances in determinant evaluation 
could facilitate sampling a very large number of projectors 

An apparent disadvantage of SHDMC is that the method 
is recursive. This disadvantage is partially removed since i) 
the number of blocks M used to collect data is increased only 
if necessary to improve the wave-function significantly ii) 
and, the propagation to large imaginary times is avoided by 
using precisely this recursive approach that accumulates the 
propagation in successive blocks. In addition, a small value 
of k St limits large fluctuations in the weights, which recently 
have been claimed to cause an exponential cost in the conver- 
gence of DMC results. 

The dominant cost of the present algorithm to obtain the 
wave-functions and their nodes scales as x rimax x nt x 
list, with n„iax being the number of excited states, ni, the 
number of projectors ^„(R) sampled, and Ust the total num- 
ber of SHDMC steps. Of course, the error and the cost depend 
on the quality of the method used to construct $„ (R) and the 
quality of the initial trial wave-functions. Systematic errors 
decrease when nf, is large, and the statistical error decreases 
when Ust increases. For a fixed absolute error, Tif, is expected 
to increase exponentially with the number of electrons 

Note that in order to describe an arbitrary wave-function of 
a system with electrons and a typical size L in D > 1 
dimensions with a resolution Rg, one needs approximately 
{L/RsY^ basis functions. The nodal surface alone re- 
quires {L/Rs)^-^ We-i) degrees of freedom. Therefore, find- 
ing an algorithm to obtain the nodes Sn (R) of any eigenstate 
n with an arbitrary interaction in a time polynomial in is 
potentially a "Philosopher's Stone" quest. However, if expo- 
nential factors actually control the accuracy of the DMC ap- 
proach, as claimed,™ just a rock solid method to find the nodes 
which simultaneously improves the wave-function (reducing 
the population fluctuations) could be considered a satisfactory 
solution. The presented work could be the basis of such a 
method. 

In ongoing work, SHDMC methods are being developed 
and tested in larger systems. 
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